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Abstract. The uniform longitudinal flow is characterized by a linear longitudinal velocity field u x (x,t) = a(t)x, where 
a(t) = ao/(l +ao0 is the strain rate, a uniform density n(t) <* a(t), and a uniform granular temperature T(t). Direct simulation 
Monte Carlo solutions of the Boltzmann equation for inelastic hard spheres are presented for three (one positive and two 
negative) representative values of the initial strain rate oq. Starting from different initial conditions, the temporal evolution 
of the reduced strain rate a* °c ao/VT, the non-Newtonian viscosity, the second and third velocity cumulants, and three 
independent marginal distribution functions has been recorded. Elimination of time in favor of the reduced strain rate a* 
shows that, after a few collisions per particle, different initial states are attracted to common "hydrodynamic" curves. Strong 
deviations from Maxwellian properties are observed from the analysis of the cumulants and the marginal distributions. 
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INTRODUCTION 

The dynamical properties of granular gases are in general much more complex than those of conventional molec- 
ular gases due to several causes such as, for instance, collisional inelasticity, frictional effects, polydispersity, non- 
sphericity, or influence of the interstitial fluid. In order to isolate the first effect from the other ones, a favorite model 
of a granular gas consists of an ensemble of identical, smooth, inelastic hard spheres with a constant coefficient of 
normal restitution a [1, 2, 3, 4]. In the dilute regime, a kinetic theory approach based on the Boltzmann equation 

«9,/(r,v,f)+vV/(r,v,f)=7[v|/,/] (1) 

has proven to be very powerful. In Eq. (1), /(r, v,f) is the one-body velocity distribution function and J[\\f,f] is the 
Boltzmann operator for inelastic collisions [5], 

In this work, we consider this simple model of a granular gas under conditions of uniform longitudinal flow 
(ULF) and analyze the temporal evolution of the velocity distribution function and its first few moments in the 
hydrodynamic stage [6], i.e., once the kinetic stage (strongly sensitive to the initial state) has decayed. The ULF 
[6, 7, 8, 9, 10, 11, 12, 13, 14] is characterized by a linear longitudinal velocity field, a uniform density, and a uniform 
granular temperature T(t): 

u x (x,t) — a(t)x, n(t) = —a(t), a(t) = - — - — , (2) 
arj 1 + a^t 

where a(t) is the strain rate. It is important to note that the constant ciq (initial strain rate) can be either positive 
(expansion of the gas) or negative (compression of the gas). The ULF is schematically depicted in Fig. 1. 
The energy balance equation is given by [6, 14] 

d t T(t) = -ja(t)T x (t)-t;(t)T{t), (3) 

where T x = P xx /n is the anisotropic temperature along the x direction (related to the normal stress P xx ) and £(t) is the 
cooling rate, which vanishes for elastic collisions (a — 1). If a(t) > (expansion), both terms on the right-hand side of 
Eq. (3) are negative and thus the granular gas monotonically loses kinetic energy, i.e., d,T(t) < 0. On the other hand, 
if a(t) < (compression) the viscous heating term ^\a(t)\T x (t) competes with the inelastic cooling term £(t)T(t) and, 
depending on the initial state, the temperature either grows or decays until a steady state is eventually reached. 
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FIGURE 1. Sketch of the ULF for (a) a(t) < and (b) a{t) > 0. 



The relevant control parameter of the problem is the reduced strain rate (which plays the role of the Knudsen 
number) 
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is an effective collision frequency, a and m being the diameter and mas of a sphere, respectively. Obviously, |a*(f)| 
increases (decreases) with time in cooling (heating) situations and reaches a stationary value a* < only if ao < 0. 

In Ref. [6] we presented direct simulation Monte Carlo (DSMC) results for the evolution of the (reduced) strain rate 
a*(t ) and the (reduced) non-Newtonian viscosity 
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for a wide ensemble of initial conditions. Parametric plots of 77* (r ) versus a* (f ) showed that, after a first (kinetic) stage 
lasting a few collisions per particle, the system reached a second (hydrodynamic) stage where, regardless of the details 
of the initial state, the curves were attracted to a common smooth "universal" curve ri*(a*). On the other hand, the 
viscosity 77* involves second-order moments only and thus the possibility that the underlying full velocity distribution 
function might still be affected by the initial preparation, even when r\*(a*) exhibits a hydrodynamic behavior, was 
not addressed in Ref. [6]. The aim of the present work is to clarify this issue by extending that analysis to higher-order 
moments, namely the second (02) and third (03) cumulants, and to the velocity distribution itself. 



UNIFORM LONGITUDINAL FLOW 

As said above, the ULF is defined by the macroscopic fields (2), together with VT = and the balance equation (3). 
At a more basic level, the velocity distribution function /(r,v,f) becomes spatially uniform when the velocities are 
referred to a Lagrangian frame moving with the flow, i.e., 

f(r,y,t)=n(t)p{V,t), V = v-u(x,f), (7) 

where p(V,f) is the probability density function and V is the peculiar velocity. After simple algebra, Eq. (1) can be 
rewritten as [6, 12] 

d T p(V,T)-ao^-[V x p(V,T)}=noJ[V\p,p], T =Hi±M. (8) 
dv x ao 

Equation (8) shows that the original ULF problem can be mapped onto the equivalent problem of a uniform gas with 
a velocity distribution «oP(V,t) and subject to the action of a non-conservative force — maoV x x. Moreover, in the 
mapped problem the temporal evolution is monitored by the scaled variable x — J dt'n(t')/no, which is unbounded 
even if qq < since in that case z — > °° when t — > |ao| _1 . In terms of p(V,f), the longitudinal temperature T x and the 
average temperature T are defined as 



UT)=m{Vl), T(t) = ^(V 2 ), ( V (V)) = |rfV V (V)p(V,T). 



(9) 



The transverse temperatures T y and T z can be denned similarly to T x . Note that T = 5 (T x + T y + T z ). The energy balance 
equation (3) can be equivalently written in the form 



a t r(T) = -|aor,(T)-Cb(T)r(T), = J d\v 2 j[y\p, P }, (io) 

where the cooling rate Co(t) in the mapped problem is related to the cooling rate £(f) of the original problem by 
£o( T ) = »oC(0/ n (0- 1 R me mapped description, the roles of strain rate and collision frequency are played by oq and 
Vo(t) = «oV(/)/n(f), respectively. Therefore, the reduced strain rate (4) remains the same in both descriptions. 

Apart from the second-order moments T x (z) and T(z), higher-order moments provide information about the 
distribution p(V,r). In particular, deviations from a Maxwellian can be characterized by the second and third 
cumulants defined as 

(V 4 ) (V 6 ) 

a 2 (r) = — — 2-1, a 3 (T) = — T + l + 3a 2 (T). (11) 

15[r(r)/m] 2 105[r(r)/m] 3 

In general, the probability distribution function p (V, t) depends on the three components of V. It is then convenient 
to introduce the marginal distributions 

/co /'Co ^00 j-oa 
dV y dV zP (\,T), Py(V y ,z)= dV x dV zP (\,z), (12) 
-00 J CO J CO J CO 

P(V,T)=V 2 JdVp(Y,T). (13) 

While the functions p x (V Xl t) and p y (V y , t) provide information about the anisotropy of the state, P(V) is the probabil- 
ity distribution function of the magnitude of the peculiar velocity, regardless of its orientation. 



UNSTEADY HYDRODYNAMIC BEHAVIOR 

As is well known, hydrodynamics is one of the key properties of normal fluids. Let us consider a gas of elastic 
particles in an arbitrary initial state defined by a certain distribution function f°(r, v). The standard evolution scenario 
starting from that initial state occurs along two consecutive stages [15]. First, during the so-called kinetic stage, 
the velocity distribution function /(r,v,f|/°), which functionally depends on the initial state, experiments a quick 
relaxation (lasting of the order of a few collisions per particle) toward a "normal" form where all the spatial and 
temporal dependence takes place through a functional dependence on the hydrodynamic fields n, u, and T, i.e., 
/(r,v,/|/°) — > f[\\n,u,T]. Next, during the hydrodynamic stage, a slower evolution occurs. While the first stage 
is very sensitive to the initial preparation of the system, the details of the initial state are practically "forgotten" in the 
hydrodynamic regime. 

An extremely important issue is whether or not the above two-stage scenario maintains its applicability in the 
inelastic case. For the sake of concreteness, let us consider the Boltzmann equation for a driven homogeneous granular 
gas (in the Lagrangian frame): 

d t f(Y,t)-(0&f(V,t)=J[V\f,f], (14) 

where & is an (isotropic or anisotropic) operator representing the external driving and CO is a constant measuring the 
strength of the driving. The operator is assumed to preserve total mass and momentum. The original problem can 
indeed be homogeneous [16, 17] or become equivalent to a homogeneous problem after a certain change of variables. 
The latter situation happens for the uniform shear flow (USF) [6, 18, 19, 20] and the ULF [see Eq. (8)]. 

Given an operator the solution to Eq. (14) depends functionally on the initial distribution f° and parametrically 
on the value of the strength co. Since the only time-dependent hydrodynamic variable is the temperature T(t), the 
existence of a hydrodynamic regime implies that, after a certain number of collisions per particle, 

/(V,f|/°,a))^n[m/2r(f)] 3/ V*(C(f);0)*(f)), C(f) = =X= , (0*(t)= (15) 
J\ \J 1 yi WJ J v v h \)h \) ^2T{t)/m K[T(t)} 7 

Here, C is the (peculiar) velocity in units of the (time-dependent) thermal speed and co* is the reduced driving 
strength, where the choices of the constant K and the exponent 7 are dictated in each case by dimensional analysis. 



The scaled velocity distribution function f*(C;(0*) should be, for a given value of the coefficient of restitution a, 
a universal function in the sense that it is independent of the initial state f° and depends on the driving strength ft) 
through the reduced quantity ft)* only. In other words, if a hydrodynamic description is possible, the different solutions 
f(y,t\f°, ft)) of the Boltzmann equation (14) would be "attracted" to the universal form (15). This has been confirmed 
by DSMC simulations for a stochastic white-noise driving (J£" = dy) at the level of the cumulants a 2 and 03 [17], for 
the USF = V y dv x ) at the level of the viscosity, the viscometric functions, and the marginal distributions [6, 20], and 
for the ULF = dy x V x ) at the level of the viscosity [6]. We conjecture that, in contrast to some claims [17], Eq. (15) 
applies as well to the case of a Gauss's driving (J? = —dy ■ V) [16] and to a combination of the latter and the stochastic 
white-noise driving [21, 22]. 

Translated to the ULF case, Eq. (15) implies that 

77* (T|p%o)^J7>*(T)), a 23 (T\ P ,a Q )^a 2 j(a*(T)), (16) 

pxy(Vxy,T|p%o) V™/2T(T)g XJ (C x , y (t); a* (t)) , P(V,T\p Q ,a Q ) -> y/m/2T{x)F{C(x);a* {%)). (17) 

As said above, the validity of the first term of Eq. (16) was addressed in Ref. [6]. In the next section we extend the 
analysis to 02(0*), a^(a*), g x (C x ;a*), g y (C y ;a*), andF(C;a*). 



RESULTS 

We have numerically solved the Boltzmann equation (8) by the DSMC method for three values of the coefficient of 
restitution (a = 0.5, 0.7, and 0.9) and three values of the strain rate (a /vo(0) = -11.26, -0.011, and 0.011). The 
two negative values of oq correspond to a compressed ULF, so that the viscous heating term | |ao | T x competes with the 
cooling term £ T in the first equation of Eq. (10). The magnitude of a /vo(0) = -11.26 is large enough as to make 
the viscous heating initially prevail over the inelastic cooling (even for a — 0.5). As the granular gas heats up, the 
cooling term grows more rapidly than the heating term until eventually both terms cancel each other and a steady state 
is reached. Conversely, the magnitude of ao/vo(0) = —0.01 1 is so small that the viscous heating is initially dominated 
by the inelastic cooling (even for a = 0.9) and the granular gas cools down. Now, the cooling term decays more rapidly 
than the heating term until the same steady state as before is eventually reached. On the other hand, the positive value 
a o/vo(0) = 0-01 1 corresponds to an expanded ULF and both terms |ao7x an d Co^ produce a cooling effect. Therefore, 
T(x) monotonically decreases (and thus a* (t) monotonically increases) without any bound and no steady state exists. 

For each one of the nine pairs (a,ao) we have considered five initial conditions. First, we have taken the local 
equilibrium state 

p°<v> = (^) 3/ V^, 

where T° is the (arbitrary) initial temperature. The initial longitudinal temperature, cumulants and marginal distribu- 
tions are simply 

7i(0) = r°, fl 2 (0)=fl 3 (0)=0, (19) 
^ = ^- mV ' ,1T \ M^^e-^, P(V,0)=4nV^^y /2 e-^ (2Q) 
Besides, we have considered four anisotropic initial conditions of the form 

P ° (V) = 2"\/2^ e ^^ (21) 
where V° = ^TT^] m is the initial thermal speed and <j> = 0, 71 / 4, n/2, and 3n/4. In this case, 

r x (0)=2r°cos 2 </>, fl2 (0) = -^, fl 3 (0) = -^, (22) 

px(V x ,0) = \[8 (V x - V°cos 0) + 8 (V x + V°cos <j>)] , p y (V y ,0) = [8 (V y - V°sin</>) + 8 (V y + V°sin</>)] , 

P(V, ) = ./3^ e -(^ 2 /2r°-i) 2V &(V 2 -2T°/m), 
V ; \ 27tT° JV 2 -2T°/m V ' ; 



(23) 
(24) 




FIGURE 2. Plot of (a) the reduced viscosity T7*, (b) the second cumulant 02, and (c) the third cumulant 03 versus the reduced 
strain rate a* for a = 0.5, a = 0.7, and a = 0.9. The circles represent the steady-state points, while the triangles represent the 
values in the HCS. 



where © is the Heaviside step function. 

In the course of the simulations we have measured T(x), T x (z), 02(^)1 aj(T), p x (V x ,T), p v (V v ,T/), and P(V,t) for 
each one of the 45 cases described above. From T(t) and T x {x) the temporal evolution of the reduced strain rate 
a*(t) [cf. Eq. (4)] and the reduced viscosity TJ*(t) [cf. Eq. (6)] has been followed. Elimination of time between both 
quantities allows one to get a parametric plot of 77* versus a*. In Ref. [6] we observed that, after a few collisions 
per particle, the curves corresponding to the five initial conditions for each one of the nine values of the pair (a,ao) 
collapse to a common hydrodynamic curve. For instance, in the case a = 0.5 the duration of the kinetic stage was 
about 3^1 collisions per particle for ao/vo(0) = —1 1-26 and about 7-8 collisions per particle for ao/vo(0) = ±0.01 1, 
while the total relaxation period toward the steady-state values a* and rj* took about 20 collisions per particle for 
flo/v (0) = -11.26 and a /v (0) = -0.011. 

Although the ULF non-Newtonian viscosity was analyzed in Ref. [6], for the sake of completeness we show in 
Fig. 2(a) rj*(a*) for three windows of a* where the hydrodynamic regime is practically established: —2 < a* < a* 
(corresponding to a /vo(0) = -11.26), -0.08 > a* > a* (corresponding to a /vo(0) = —0.011), and a* > 0.08 
(corresponding to ao/Vo(0) = 0.011). A non-monotonic behavior of rj*(a*) is observed, with a maximum at a* = 
—0.47 (a = 0.5 and 0.7) and a* — —0.44 (a — 0.9). In the figure the open circles represent the steady state points 
(a* , 77 *) and the open triangles at a* = represent values of the Navier-Stokes (NS) viscosity obtained independently 
[23, 24, 25, 26, 27]. While, for each a, the steady-state point is an attractor of the heating (a* < a*) and cooling 
(a* > a*) branches with negative a*, the NS point is a repeller of the two cooling branches (a* < a* < and a* > 0) 
[13, 14]. It is noteworthy that, although the need of discarding the kinetic stage creates the gap —0.08 > a* > 0.08, the 
extrapolations of the cooling branches with positive and negative a* smoothly join at the NS point. 

In the case of the heating states (ao/vo(0) = —11 .26) the evolution is so fast when starting from the local equilibrium 
initial condition (18) and from the anisotropic initial condition (21) with = 7r/2 that the corresponding curves join 
the hydrodynamic line past the maximum located at a* < a* . Thus, in Fig. 2 and henceforth those two initial conditions 
are omitted in the case ao/ Vo (0) = — 1 1 .26. 

Figures 2(b) and 2(c) show the cumulants 02(0*) and 03(0*), respectively, for the same windows of a* as in Fig. 
2(a). Again, the open circles represent the steady-state points. The open triangles at a* = correspond to homogeneous 
cooling state (HCS) values obtained independently [16, 28, 29, 30]. We observe that, for each a, the HCS values are 
fully consistent with the extrapolation to a* — >• of the two cooling branches. For these higher-order moments, the 
overlapping of the individual heating curves in the region —2 < a* < — 1 is less robust than in the case of 77*. As 
indicated before, this is a consequence of the very fast evolution taking place for the states with ao/vo(0) = —11.26. 
In fact, the value a* = —2 is reached after less than 0.3 collisions per particle only. We have checked that starting 
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FIGURE 3. Marginal probability distribution functions (a) g x (C x ), (a)g y (C y ), and(c)F(C) for a = 0.5 at a* = -0.5 and a* =0.5. 
The thin solid lines represent the corresponding Maxwellian distributions 
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FIGURE 4. Marginal probability distribution function F(C) for a = 0.5, 0.7, and 0.9 at (a) a* = -0.5 and (b) a* = 0.5. 



from more negative values of ao/vo(0) hardly changes the situation. We also observe that, for each a, the noise in the 
cooling curves with positive a* increases as the order of the moment grows. 

The dependence of «2 on a* is similar for the three values of a. The second cumulant is negative at a* = —2, grows 
with increasing a*, changes sign, and reaches a maximum value at a* — —0.45 (a = 0.5), a* = —0.40 (a = 0.7), 
and a* — —0.33 (a = 0.9). Thereafter, ai decays, reaches a local minimum at a* = and then grows with increasing 
positive a*. At a fixed value of a*, we observe that «2 increases with increasing inelasticity. The dependence of 123 
on a* for compression states (a* < 0) is more complex than that of 02- Instead of a maximum, —03(0*) presents a 
minimum at a* = -0.70 (a = 0.5), a* = -0.61 (a = 0.7), and a* = -0.48 (a = 0.9), followed by a maximum at 
a* = —0.41 (a = 0.5), a* = —0.36 (a = 0.7), and a* = —0.25 (a = 0.9). Moreover, the curves corresponding to the 
different values of a cross each other in the region of negative a*. 

Now we turn our attention to the marginal distributions g x (C x ;a*), g y (C y ;a*), and F{C;a*) [cf. Eq. (17)]. Although 



we have evaluated those functions for the whole temporal evolution of the systems, here we fix two representative 
values of the reduced strain rate: a* = —0.5 and a* = 0.5. The first value belongs to the heating compression branch 
and is reached after slightly less than 2 collisions per particle. The second value belongs to the cooling expansion 
branch and is reached after 13 (a — 0.5), 18 (a — 0.7), or 38 (a = 0.9) collisions per particle. The functions g x 
and g y for a — 0.5 are plotted in Figs. 3(a) and 3(b), respectively. The representation in the horizontal and vertical 
axes is chosen to visualize deviations from the Maxwellians g%y(C Xt y) — (7zT xy /T)~ l / 2 exp(—C xy T/T xy ). We first 
note a high degree of collapse of data corresponding to different initial conditions to a common curve in the velocity 
domain C xy < 6. Comparison of g x (C x ) for a* — —0.5 and a* = 0.5 shows a much broader distribution in the first case 
than in the second one. The opposite happens for g y (C y ). This is consistent with the measured values T x /T = 2.09 and 
Ty/T = 0.45 at a* = -0.5 and T x /T = 0.49 and T y /T = 1 .26 at a* = 0.5. Apart from that, we observe an overpopulated 
high-velocity tail with respect to the Maxwellian gf y (C xy ). In the case of the broader distributions this overpopulation 
phenomenon occurs beyond the region C xy < 6 displayed in Figs. 3(a) and 3(b). 

The probability distribution function for the reduced speed, F(C), is plotted in Fig. 3(c) for a = 0.5. Again, the 
representation is chosen to highlight deviations from the Maxwellian F M (C) — An~ l l 2 C 2 e~ c . We see that F(C) for 
both the compression (a* = —0.5) and the expansion (a* = 0.5) states exhibits strong departures from the Maxwellian 
F M (C), with overpopulated tails. This is especially so in the case of a* = —0.5, in agreement with the fact that and 
|«3 1 are clearly larger at a* = —0.5 than at at a* — 0.5 [cf. Figs. 2(b) and 2(c)]. 

The function F(C) is plotted in Fig. 4 for the three values of the coefficient of restitution and the two chosen 
values of the reduced strain rate. The main observation is that F(C) for a* = —0.5 presents a singular behavior at 
C 2 s=s 12, 7.5, and 5 for a = 0.5, 0.7 and 0.9, respectively. This is a remaining artifact associated with the special 
initial conditions (21). As shown by Eq. (24), at t — all the particles have a speed V > y/2T°/m and the distribution 
function diverges at V > y/2T°/m. Since a* = —0.5 is reached after not more than 2 collisions per particle, there 
exists a certain population of surviving particles which have not collided yet. Those particles are subject to the heating 
effect due to the non-conservative external force —maoV x x but not to the collisional cooling effect. Consequently, they 
increase their energy more than the rest of the particles, have an ever increasing reduced speed C and thus contribute to 
the tail of F (C) only. As the inelasticity decreases the cooling effect becomes less important and the reduced speed of 
the surviving particles grows more slowly. Therefore, in the case a* = —0.5, the tail of the distribution function F(C) 
for values of C equal to or larger than the singularity is still dependent on the initial state and cannot be considered as 
hydrodynamic yet. On the other hand, the tail does not have any practical influence on the first few velocity moments. 
In the case a* = 0.5 the number of collisions is much higher and thus the influence of the surviving particles is 
negligible. 

CONCLUSIONS 

In this paper we have presented DSMC numerical solutions of the inelastic Boltzmann equation for the unsteady ULF 
in the co-moving Lagrangian frame. In order to uncover the three types of possible regimes (heating compression 
states, cooling compression states, and cooling expansion states), three values of the initial strain rate have been 
considered. Starting from different initial conditions, the temporal evolution of the granular temperature T, the 
longitudinal temperature T x , the velocity cumulants ai and aj, and the marginal probability distribution functions 
Px(V x ), py(V y ), and P(V) has been recorded. By eliminating time in favor of the reduced strain rate a* we have checked 
that, after a first kinetic stage lasting a few collisions per particle, the curves corresponding to different initial states 
tend to collapse to common "hydrodynamic" curves. These findings extend to the cumulants and to the distribution 
function recent results [6] obtained for the non-Newtonian viscosity. The dependence of the cumulants ai and 03 
on a* exhibits a non-monotonic behavior with interesting features. In consistency with this, the velocity distribution 
functions are seen to depart from a Maxwellian form. 
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